---
jupytext:
  formats: md:myst,ipynb
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.4
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

```{code-cell}
:init_cell: true

import mmf_setup

mmf_setup.nbinit()
```

# Travelling Waves

+++

A traveling wave with velocity $v$ has the following form:

$$
  \psi(x, t) = \sqrt{\rho}(x - vt) e^{\I\phi(x-vt) - \I\omega t}, \\
  \I\hbar \dot{\psi} = -\frac{\hbar^2\psi''}{2m} + (g\rho - \mu) \psi.\\
  -\I\hbar v\psi' + \hbar\omega\psi = -\frac{\hbar^2\psi''}{2m} + (g\rho - \mu) \psi.
$$

By appropriately adjusting $\mu$ we may set $\omega = 0$, hence we may drop the factor $e^{-\I\omega t}$ from the wavefunction:

$$
  \psi(x - vt) = \sqrt{\rho}(x - vt) e^{\I\phi(x-vt)}.
$$

The last equation above now represents a stationary solution in a moving frame with coordinate $y = x-vt$:

$$
  0 = -\frac{\hbar^2\psi''(y)}{2m} + \I\hbar v\psi'(y) + [g\rho(y)- \mu]\psi(y)
             = \left[\frac{\op{p}^2}{2m} -\op{p} v + g\rho(y) - \mu\right]\psi(y)
             = \left[\frac{(\op{p} - p_v)^2}{2m} - \frac{p_v^2}{2m} + g\rho(y) - \mu\right]\psi(y)
$$

where $p_v = mv$.  Finally, we may recast this in terms of the original GPE by transforming $\psi(y)$ as follows:

$$
  \psi(y) = e^{\I p_v x/\hbar}\tilde{\psi}(y)\\
  (\op{p}-p_v)\psi(y) = e^{-\I p_v x/\hbar}\op{p}\tilde{\psi}(y)\\
  \left[\frac{\op{p}^2}{2m} + g\rho(y) - \left(\mu + \frac{p_v^2}{2m}\right)\right]\tilde{\psi}(y)
$$

If the original function $\psi(y+L) = \psi(y)$ was periodic, then $\tilde{\psi}(y+L) = e^{\I p_v L/\hbar}\tilde{\psi}(y)$.  In other words, twisted boundary conditions should be applied with a twist $\theta = p_vL/\hbar$.

+++

# GPE

+++

$$
  \psi(x-vt) = e^{\I\phi(x-vt)}f(x-vt), \qquad f(x) = \sqrt{\rho(x)}\\
  \psi' = \I\phi'\psi + e^{\I\phi}f', \qquad
  \psi'' = \I\phi''\psi - (\phi')^2\psi +2\I\phi'\psi'
  + e^{\I\phi}\I\phi'f' + e^{\I\phi}f'',\\
  \I\dot{\psi} = -v\I\psi' = v\phi'\psi - v\I e^{\I\phi}f'
$$
$$
  2v\phi'\psi - 2v\I e^{\I\phi}f' =
  \I\phi''\psi - (\phi')^2\psi +2\I\phi'\psi' + e^{\I\phi}\I\phi'f' + e^{\I\phi}f'' + 2e^{\I\phi}f^3\\
  [3\phi'+2v]f' = -\phi''f\\
  [3\phi' + 2v]\phi'f = f'' + 2f^3\\
  -\left(\frac{(\phi')^2}{2}\right)'
  = \frac{f'f''}{f^2} + (f^2)'
  = \frac{f'f''}{f^2} + \rho'\\
$$

$$
  f^2 = \rho, \qquad
  2ff' = \rho', \qquad
  2f'f' + 2ff'' = \rho'', \qquad
$$

+++

Here we consider the problem of finding traveling wave solutions in a BEC.  From a numerical perspective, it is highly beneficial if the problem can be stated in terms of a well-defined minimization problem.  To start, we consider the available analytic solution for the conventional GPE.  These solutions are presented in [El:2016]:

$$
  \DeclareMathOperator{\sn}{sn}
  \I\dot{\psi} = \frac{-\psi''}{2} + \abs{\psi}^2\psi\\
  \rho = \abs{\psi}^2 = \frac{(r_4-r_3-r_2+r_1)^2}{4}
  + (r_4-r_3)(r_2-r_1)\sn^2(\sqrt{(r_4-r_2)(r_3-r_1)}\xi; m)\\
  C = \frac{(-r_1-r_2+r_3+r_4)(-r_1+r_2-r_3+r_4)(r_1-r_2-r_3+r_4)}{8}\\
  u = v - \frac{C}{\rho}, \qquad
  \xi = x-vt - \xi_0, \qquad
  v = \frac{1}{2}\sum r_i, \\
  a = (r_4-r_3)(r_2-r_1), \qquad
  r_1 \leq r_2 \leq r_3 \leq r_4, \qquad
  m = \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)},\\
  L = \frac{1}{2}\oint \frac{\d\lambda}
    {\sqrt{(\lambda - r_1)(\lambda - r_2)(\lambda - r_3)(\lambda - r_4)}}
    = \frac{2K(m)}{\sqrt{(r_4-r_2)(r_3-r_1)}}.
$$

The physical interpretations are:

* $a$: Amplitude
* $v$: Phase velocity
* $u$:

+++

$$
  \DeclareMathOperator{\sn}{sn}
  \I\dot{\psi} = \frac{-\psi''}{2} + \abs{\psi}^2\psi\\
  \rho = \abs{\psi}^2 = \alpha + a\sn^2(z; m)\\
  \alpha = \frac{(r_4-r_3-r_2+r_1)^2}{4}, \\
  D = \sqrt{(r_4-r_2)(r_3-r_1)}, \qquad
  z = D\xi\\
  C = \frac{(-r_1-r_2+r_3+r_4)(-r_1+r_2-r_3+r_4)(r_1-r_2-r_3+r_4)}{8}\\
  u = v - \frac{C}{\rho}, \qquad
  \xi = x-vt - \xi_0, \qquad
  v = \frac{1}{2}\sum r_i, \\
  a = (r_4-r_3)(r_2-r_1), \qquad
  r_1 \leq r_2 \leq r_3 \leq r_4, \qquad
  m = \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)},\\
  L = \frac{1}{2}\oint \frac{\d\lambda}
    {\sqrt{(\lambda - r_1)(\lambda - r_2)(\lambda - r_3)(\lambda - r_4)}}
    = \frac{2K(m)}{\sqrt{(r_4-r_2)(r_3-r_1)}}.
$$

+++

$$
  \DeclareMathOperator{\cn}{cn}
  \DeclareMathOperator{\dn}{dn}
  \sn'(z) = \cn(z)\dn(z), \qquad
  \cn'(z) = -\sn(z)\dn(z), \qquad
  \dn'(z) = -k^2\sn(z)\cn(z)\\
  \sn''(z) = -\sn(z)[\dn^2(z) +k^2\cn^2(z)], \qquad
  \cn'(z) = -\sn(z)\dn(z), \qquad
  \dn'(z) = -k^2\sn(z)\cn(z)\\
  \cn^2 + \sn^2 = \dn^2 + k^2\sn^2 = 1.
$$

A special limit is when $m=1$:

$$
  \sn(z;m=1) = \tanh(z), \qquad
  \cn(z;m=1) = \dn(z;m=1) = \sech(z).
$$

+++

$$
  \psi = \sqrt{1+\alpha\sn^2(z)}\\
  \psi' = \frac{\alpha\sn(z)\sn'(z)}{\psi}\\
  \psi'' = \frac{\alpha}{\psi}\left(
  \sn'(z)\sn'(z) + \sn(z)\sn''(z)
  -\frac{\alpha\sn^2(z)[\sn'(z)]^2}{\psi^2}\right)\\
$$

+++

$$
  \psi' = \frac{\alpha\sn(z)\cn(z)\dn(z)}{\psi}\\
  \psi'' = \frac{\alpha}{\psi}\left(
  \cn^2(z)\dn^2(z) - \sn^2(z)\dn^2(z) - k^2\sn^2(z)\cn^2(z)
  -\frac{\alpha\sn^2(z)\cn^2(z)\dn^2(z)}{\psi^2}\right)\\
$$

+++

First we test these.  To match units we set $\hbar = m = g = 1$.

```{code-cell}
K(0.9)
```

```{code-cell}
:init_cell: true

from gpe.imports import *
from scipy.special import ellipj, ellipk
import gpe.bec


def sn(u, m):
    return ellipj(u, m)[0]


def K(m):
    return ellipk(m)


class State(gpe.bec.State):
    def __init__(self, Nx=32, rs=[0.1, 0.2, 0.3, 0.4], xi0=0):
        r1, r2, r3, r4 = rs = sorted(rs)
        k = np.sqrt((r4 - r2) * (r3 - r1))
        m = (r2 - r1) * (r4 - r3) / (r4 - r2) / (r3 - r1)
        v = sum(rs) / 2.0
        L = 2.0 * K(m) / k
        gpe.bec.State.__init__(self, Nxyz=(Nx,), Lxyz=(L,), m=1.0, g=1.0, hbar=1.0)
        x = self.xyz[0]
        self.T = abs(L / v)
        t = 0
        xi = x - v * t - xi0
        a = (r4 - r3) * (r2 - r1)
        # a = m**2*k**2
        # rho_0 = (r4-r3-r2+r1)**2/4.0
        rho_0 = (r4 - r3 + r2 - r1) ** 2 / 4.0
        rho = rho_0 + a * sn(k * xi, m) ** 2
        self.a_ = self.a = a
        self.m_ = m
        self.k_ = self.k = k
        self.rho_0_ = self.rho_0 = rho_0

        k_ = self.m * v / self.hbar
        self[...] = np.exp(1j * k_ * x) * np.sqrt(rho)

    def get_Vext(self):
        return 0.0


s = State(Nx=256, rs=[-2.0, -1.0, 1.0, 2.0])
```

```{code-cell}
n = s.get_density()
x = s.xyz[0]
x_ = x[1:-1]
psi_ = s[1:-1]
ddpsi_ = np.diff(np.diff(s[...])) / np.diff(x)[1:] ** 2
plt.plot(x_, ddpsi_ / psi_ / 2 + abs(psi_) ** 2)
```

Check the soliton limit (2.122)

```{code-cell}
s = State(Nx=128, rs=[-4.000001, 1.0, 1.000001, 2.0])
x = s.xyz[0]
s.plot()
rho_min = abs(s[...]).min() ** 2
rho_max = abs(s[...]).max() ** 2
plt.axhspan(rho_min, rho_min + s.a, fc="y", alpha=0.5)
# plt.plot(x, rho_max-s.a/np.cosh(np.sqrt(s.a)*x)**2, '+:')
```

```{code-cell}
rho_0 = s.rho_0
s[...] = np.sqrt(rho_0) * np.tanh(np.sqrt(rho_0) * x)
s.plot()
```

```{code-cell}
s.cooling_phase = 1
e = EvolverABM(s, dt=0.2 * s.t_scale)
s = State(Nx=256, rs=[-2.0, -1.0, 1.0, 2.0])

with NoInterrupt(ignore=True) as interrupted:
    while e.t < e.y.T and not interrupted:
        e.evolve(1000)
        plt.clf()
        e.y.plot()
        display(plt.gcf())
        clear_output(wait=True)
```

```{code-cell}
%debug
```

# Solitons

+++

The GPE admits grey solitons with infinite period.  These solitons arise from the previous solutions in the limit where $m\rightarrow 1$ so that $\sn\rightarrow\tanh$ and $\cn\rightarrow\dn\rightarrow \sech$.

The soliton solution is:

$$
  \psi(x, t) = \sqrt{\bar{\rho}}e^{-\I\mu t/\hbar}\left[
  \I\frac{v}{c}
  + \frac{u}{c}\tanh\left(\frac{x-vt}{l}\right)\right], \qquad
  \rho = \bar{\rho}\left[\frac{v^2}{c^2}
  + \frac{u^2}{c^2}\tanh^2\left(\frac{x-vt}{l}\right)\right]\\
  u^2+v^2 = c^2, \qquad \mu=g\bar{\rho} = mc^2,\qquad l=\hbar/m u.
$$

This can be expressed as

$$
  \psi(x, t) = \sqrt{\rho}e^{-\I\mu t/\hbar + \I\phi(x)}, \qquad
  \tan\phi(x) = \frac{v}{u\tanh\left(\frac{x-vt}{l}\right)}, \qquad
  \phi_{\mathrm{twist}} = \phi(\infty) - \phi(-\infty) = 2\tan^{-1}\frac{v}{u}
$$

+++

$$
  \rho(0) = \bar{\rho}\frac{v^2}{c^2}
$$

+++

or, in the same units $m=\hbar = g = 1$ as El and Hoefer:

$$
  \rho = \bar{\rho}\left[\frac{v^2}{c^2} + \left(1-\frac{v^2}{c^2}\right)\tanh^2[u(x-vt]\right]\\
  u= \sqrt{c^2-v^2}, \qquad \mu=\bar{\rho} = c^2
$$


To compare with (2.117):

$$
  \rho = \frac{(r_4-r_3-r_2+r_1)^2}{4} +(r_4-r_3)(r_2-r_1)\sn^2[\sqrt{(r_4-r_2)(r_3-r_1)}(x-Vt), m]
$$

we have $m=1$, $\sn = \tanh$ and

$$
  \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)} = m = 1,\\
  \frac{(r_4-r_3-r_2+r_1)^2}{4} = \bar{\rho}\frac{v^2}{c^2} = v^2\\
  a = (r_4-r_3)(r_2-r_1) = \bar{\rho}\left(1-\frac{v^2}{c^2}\right) = \bar{\rho} - v^2\qquad
  \sqrt{(r_4-r_2)(r_3-r_1)} = u, \\
  V = \frac{r_1+r_2+r_3+r_4}{2}
$$

+++

Everything here is reasonable except the definition of the phase velocity $V$ which differs from the soliton velocity $v$.

+++

The background density is

$$
  \bar{\rho} = \frac{(r_4-r_3-r_2+r_1)^2}{4} + (r_4-r_3)(r_2-r_1)
             = \frac{(r_4-r_3+r_2-r_1)^2}{4}.
 $$

+++

To compare with (2.122):


, the stationary solution must have $\bar{\rho} = a_s$ or:

$$
  (r_4-r_1)^2 = 4(r_4-r_2)(r_2-r_1)
$$

If $r_1 = -r_4-2r_2$, then

$$
  4(r_4+r_2)^2 = 4(r_4-r_2)(r_4+3r_2)
$$

+++

# Issues

+++

I was having some issues finding solutions so I considered stationary solutions:

$$
  \psi = \sqrt{\rho_0 + a\sn^2(kx,m)}.
$$

Symbolically solving for solutions to the GPE (in Maple) gives the following solutions:

\begin{align}
  a&=k^2m^2, &&& \mu &= \rho_0 - \frac{k^4m^2}{2\rho_0}, \\
  a&=0, &&& \mu &= \rho_0, & \psi &= \sqrt{\rho_0}.\tag{Constant}\\
  a&=\rho_0 m^2 , & k&=\sqrt{\rho_0}, & \mu &= \rho_0- \frac{\rho_0m^2}{2}, & \psi &= \sqrt{\rho_0}\sqrt{1+m^2\sn^2(\sqrt{\rho_0}x,m)}.\tag{Constant}\\
  a&=\rho_0 , & k&=\frac{\sqrt{\rho_0}}{m}, & \mu &= \rho_0 - \frac{\rho_0}{2m^2}, & \psi &= \sqrt{\rho_0}\sqrt{1+\sn^2\left(\frac{\sqrt{\rho_0}}{m}x,m\right)}.\tag{Constant}\\
\end{align}

(assuming $a,k\neq 0$ and real $k$) gives:

$$
  m^2 = \frac{a}{k^2}, \qquad
  \mu = \rho_0 - \frac{ak^2}{2\rho_0}, \qquad
  \rho = \rho_0[1 - \sn^2(kx,m)].
$$

```maple
restart;
X_:=JacobiSN(k*x,sqrt(m2));
psi:=sqrt(rho[0] + a*X_^2);
Rho:=psi^2:
res:=collect(numer(simplify(-diff(psi, x, x)/2/psi + Rho - mu)), X_):
mu:=expand(solve(coeff(res, X_, 0), mu));
m2:=solve(coeff(res, X_^6), m2);
solve([coeffs(res, X_)], [a,k^2]);
```

+++

The dark soliton solution has $m=1$, $a=k^2$

```{code-cell}
s.m
```

$$
  a = (r_4-r_3)(r_2-r_1)\\
  a = mk^2\\
  k^2 = (r_4-r_2)(r_3-r_1)\\
$$

```{code-cell}
s = State(rs=[-4.1, 1.0, 1.1, 2.0])
x = s.xyz[0]
m = s.m_
a = s.a_
k = s.k_
rho_0 = s.rho_0_
print (m, a, m ** 2 * k ** 2, rho_0, k)
k = np.sqrt(rho_0)
plt.plot(x, 1 + sn(k * x, m) ** 2)
# s[...] = np.sqrt(rho_0*(1-sn(x,m)))
```

$$
  m^2 = \frac{a\alpha^2}{D^2}, \qquad
  (3m^2 - a) = 6\frac{a\alpha^2}{D^2}, \qquad
  -2a(m^2 + 1) = 6\frac{a\alpha^2}{D^2}, \qquad
  \mu = \frac{-D^2a}{2} + \alpha^2
$$

$$
  -2a(m^2 + 1) = (3m^2 - a) = 6m^2\\
  a = -3m^2\\
  m^2 + 1 = 1
$$

+++

To do this, we first consider boosting to a moving frame with velocity $v$ so that in this frame, the traveling wave solution is stationary and periodic.  We second
Here we are looking for solutions that satisfy:

$$
  \psi(x + L) = e^{\I m v L / \hbar}\psi(x)
$$

+++

# Minimization

+++

Here we consider the following hypothesis:

The traveling waves can be found as minimum energy solutions in a box of period $L$ with twisted boundary conditions holding both the total particle number fixed and the value of the wavefunction at one point.

+++

We start with the soliton solution:

$$
  \rho = \bar{\rho}\left[\frac{v^2}{c^2} + \left(1-\frac{v^2}{c^2}\right)\tanh^2[u(x-vt]\right]\\
  u= \sqrt{c^2-v^2}, \qquad \mu=\bar{\rho} = c^2.
$$

At the core, the density is:

$$
  \bar{\rho}\frac{v^2}{c^2}
$$

```{code-cell}
from gpe.imports import *
from scipy.special import ellipj, ellipk
import gpe.bec, gpe.minimize

reload(gpe.bec)


class State(gpe.bec.StateTwist):
    def __init__(
        self,
        Nx=32,
        L=10.0,
        mu=1.0,
        psi_0=1.0,
        ind=None,
        v=0.0,
        twist=None,
        m=1.0,
        hbar=1.0,
        g=1.0,
    ):
        if ind is None:
            ind = Nx // 2
        self.ind = ind
        self.psi_0 = psi_0
        self.p_v = p_v = m * v
        if twist is None:
            twist = np.pi + p_v * L / hbar
        gpe.bec.StateTwist.__init__(
            self, Nxyz=(Nx,), Lxyz=(L,), m=m, hbar=hbar, g=g, mu=mu, twist=(twist,)
        )
        self[self.ind] = psi_0

    def exact_psi(self):
        """Return the analytic soliton"""
        v = abs(self.psi_0)
        u = np.sqrt(self.mu - v ** 2)
        return 1j * v + u * np.tanh(u * self.xyz[0])

    def get_Vext(self):
        return -(self.mu + self.p_v ** 2 / 2.0 / self.m)

    def _compute_dy_dt(self, dy, subtract_mu=False):
        dy = gpe.bec.StateTwist.compute_dy_dt(self, dy=dy, subtract_mu=subtract_mu)
        #    dy[self.ind] = 0
        return dy


class MinimizeState(gpe.minimize.MinimizeState):
    def __init__(self, state, **kw):
        self.psi_0 = state.psi_0
        self.ind = state.Nxyz[0] // 2
        gpe.minimize.MinimizeState.__init__(self, state, **kw)

    def pack(self, psi):
        psi[self.ind] = self.psi_0
        fact = 1 if self.real else 2
        x = gpe.minimize.MinimizeState.pack(self, psi)
        return np.concatenate([x[: self.ind * fact], x[(self.ind + 1) * fact :]])

    def unpack(self, x, state=None):
        fact = 1 if self.real else 2
        x = np.concatenate([x[: self.ind * fact], [0, 0], x[self.ind * fact :]])
        state = gpe.minimize.MinimizeState.unpack(self, x, state=state)
        assert np.allclose(state[self.ind], 0)
        state[self.ind] = self.psi_0
        # plt.clf()
        # state.plot()
        # plt.twinx()
        # plt.plot(state.xyz[0], np.angle(state[...]/state.twist_phase), 'r--')
        # display(plt.gcf())
        # clear_output(wait=True)
        return state
```

```{code-cell}
s = State(Nx=128 * 8, L=120.0, mu=1.0, psi_0=0.5, v=v)
psi_s = s.exact_psi()
v = (
    (np.angle(psi_s[-1]) - np.angle(psi_s[0]) + np.pi)
    * s.hbar
    / s.m
    / (s.Lxyz[0] - s.Lxyz[0] / s.Nxyz[0])
)
s = State(Nx=128 * 8, L=120.0, mu=1.0, psi_0=0.5, v=v)
s[...] = s.exact_psi()
print (v)
```

```{code-cell}
plt.plot(s[...] / s.twist_phase)
```

```{code-cell}
v
```

```{code-cell}
s = s1
s.cooling_phase = 1.0
e = EvolverABM(s, dt=0.5 * s.t_scale)

with NoInterrupt(ignore=True) as interrupted:
    while not interrupted:
        e.evolve(100)
        plt.clf()
        e.y.plot()
        display(plt.gcf())
        clear_output(wait=True)
```

```{code-cell}
s = State(Nx=128 * 4, L=20.0, mu=1.0, psi_0=0.1, v=0.1)
s[...] = 1.0
m = MinimizeState(s, fix_N=False)
s1 = m.minimize(use_scipy=True)
plt.plot(s1.xyz[0], s1.get_density() - abs(s1.exact_psi()) ** 2)
# s1.plot()
abs(s1.get_density() - abs(s1.exact_psi()) ** 2).max()
```

```{code-cell}
plt.plot(s1.xyz[0], np.angle(s1[...]))
plt.plot(s1.xyz[0], np.angle(s1.exact_psi()))
```

```{code-cell}
vs = np.linspace(0, 0.4, 10)
errs = []
for v in vs:
    s = State(Nx=128 * 4, L=20.0, mu=1.0, psi_0=0.2, v=v)
    m = MinimizeState(s, fix_N=False)
    s1 = m.minimize(use_scipy=True)
    errs.append(abs(s1.get_density() - abs(s1.exact_psi()) ** 2).max())
```

```{code-cell}
plt.plot(vs, errs, "-+")
```

```{code-cell}
s = gpe.bec.State(Nxyz=(64,), Lxyz=(10.0,))
s[...] = 1.0
m = gpe.minimize.MinimizeState(
    fix_N=False,
)
s1 = m.minimize(use_scipy=True)
s1.plot()
```

```{code-cell}
Es = []
m = MinimizeState(s)
m.check()
s1 = m.minimize(use_scipy=True, fix_N=False, callback=callback)
# plt.plot(s1.xyz[0], np.log10(abs(abs(s1[...])**2-s.mu)))
s1.plot()
plt.plot(s1.xyz[0], abs(s1.exact_psi()) ** 2)

plt.plot(Es)
```

```{code-cell}
s1.get_density().max(), 1 + (s1.k_B[0] ** 2) / 2
```

```{code-cell}
-s1.mu, s1.get_Vext()
```

```{code-cell}
Es = []
twists = np.linspace(0, 2 * np.pi, 10)
for twist in twists:
    s = State(Nx=128, L=50.0, mu=1.0, psi_0=0.05, twist=twist)
    m = MinimizeState(s)
    s1 = m.minimize(use_scipy=True, fix_N=False)
    Es.append(s1.get_energy())
    print twist, Es[-1]
```

```{code-cell}
s1.plot()
```

```{code-cell}
plt.plot(twists, Es)
```

```{code-cell}
m = MinimizeState(s)
s1 = m.minimize(use_scipy=True)
plt.clf()
s1.plot()
s1.get_energy()
```